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One of the longstanding issues in numerical relativity is to enable a simulation taking 
account of microphysical processes (e.g., weak interactions and neutrino cooling). We de- 
velop an approximate and explicit scheme in the fully general relativistic framework as a 
first implementation of the microphysics toward a more realistic and sophisticated model- 
ing. In this paper, we describe in detail a method for implementation of a realistic equation 
of state, the electron capture and the neutrino cooling in a multidimensional, fully general 
relativistic code. The procedure is based on the so-called neutrino leakage scheme. To check 
the validity of the code, we perform a two dimensional (2D) simulation of spherical stel- 
lar core collapse. Until the convective activities set in, our results approximately agree, or 
at least are consistent, with those in the previous so-called state-of-the-art simulations. In 
particular, the radial profiles of thermodynamical quantities and the time evolution of the 
neutrino luminosities agree quantitatively. The convection is driven by negative gradients of 
the entropy per baryon and the electron fraction as in the previous 2D Newtonian simula- 
tions. We clarify which gradient is more responsible for the convection. Gravitational waves 
from the convection are also calculated. We find that the characteristic frequencies of the 
gravitational-wave spectra are distributed for higher frequencies than those in Newtonian 
simulations due to the general relativistic effects. 



§1. Introduction 



1.1. Motivation 



Gravitational collapse of massive stellar core to a neutron star or a black hole 
and the associated supernova explosion are one of the important and interesting 
events in the universe. Prom observational view point, they are among the most 
energetic events in astrophysics, producing a wide variety of observable signatures, 
namely, electromagnetic radiation, neutrinos, and gravitational radiation. 

Most of the energy liberated in the collapse is eventually carried away by neu- 
trinos from the system. The total energy of neutrinos emitted is ~ GM^ s /i?NS ~ 
O.IMnsc 2 ~ several times 10 53 ergs, where Mns and -Rns are the mass and radius 
of the neutron star. Observations of gravitational collapse by neutrino detectors 
will provide important information of the deep inside of the core, because neutrinos 
can propagate from the central regions of the stellar core almost freely due to their 
small cross-sections with matters. Electromagnetic radiation, by contrast, interacts 
strongly with matters and thus gives information of collapse coming only from lower- 
density regions near the surface of the star. Bursts of neutrinos were first detected 
simultaneously by the Kamiokande^ 1 and Irvine- Michigan-BrookhaverPJ facilities in 
the supernova SN1987A, which occurred on February 23, 1987 in the Large Magel- 
lanic Cloud (for a review, see Ref. SJ)). Future detection of neutrinos will provide a 
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direct clue to reveal the physical ingredient for the supernova explosion mechanism. 

Gravitational wave astronomy will start in this decade. The first generation of 
ground-based interferometric detectors (LIGO, 4 ^ VIRGO© GEO600,® TAMA30CP) 
are now in the scientific search for gravitational waves. Stellar core collapse is one 
of the important sources for these observatories. Observations of gravitational col- 
lapse by gravitational-wave detectors will provide unique information, complemen- 
tary to that derived from electromagnetic and neutrino detectors, because gravita- 
tional waves can propagate from the innermost regions of a progenitor star to the 
detectors without attenuation by matters. Combination of the signatures of neutri- 
nos and gravitational waves will provide much information about processes of the 
core collapse and ultimately, the physics that governs the stellar core collapse. 

To obtain physically valuable information from these observations, it is neces- 
sary to connect the observed data and the physics behind it. For this purpose, a 
numerical simulation is the unique approach. However, simulating the stellar core 
collapse is one of the challenging problems because a rich diversity of physics has 
to be taken into account. All four known forces in nature are involved and play 
important roles during the collapse. General relativistic gravity plays an essential 
role in formation of a black hole and a neutron star. The weak interactions govern 
energy and lepton-number losses of the system. In particular, neutrinos transport 
most of the energy released during the collapse to the outside of the system. The 
electromagnetic and strong interactions determine the ingredient of the collapsing 
matter and the thermodynamical properties of the dense matter. Strong magnetic 
fields, if they are present, would modify the dynamics of the collapse, subsequent 
supernova explosion, and evolution of proto-neutron stars. 

Due to these complexities, the explosion mechanism of core-collapse supernovae 
has not been fully understood in spite of the elaborate effort in the past about 40 
years MM Recent numerical studiesP OOflUflSJE) have clarified that on the 
assumption of the spherical symmetry, the explosion does not succeed for the iron 
core collapse with the currently most elaborate input physics (neutrino interactions, 
neutrino transfer, and equation of states of the dense matter) on the basis of the 
standard "neutrino heating mechanism"^ (but see Ref. [T7|) for successful explosion 
in O-Ne-Mg core collapse). To increase the neutrino- heating efficiency, a wide va- 
riety of multi-dimensional effects have been explored (for recent reviews, see e.g., 
Refs. |9|), E]) and also Refs. [T8|) and [T9|) for simulations where successful explosions 
are obtained). However, it has not been completely clarified yet whether the increase 
of the heating efficiency due to such multi-dimensional effects suffices for yielding 
successful explosion, because the explosion energy resulting from these works is too 
low ~ 10 50 ergs. 

Similarly, accurate predictions of gravitational waveforms are still hampered by 
the facts that reliable estimates of waveforms require a general relativistic treat- 
ment,^ and that appropriate treatments of microphysics such as nuclear equation 
of state (EOS), the electron capture, and neutrino emissions and transfers. Indeed, 
previous estimates of waveforms have relied either on Newtonian simulations with 
including microphysics to some extent JIDJIUJIDJUIJIDJIO -ED J2EJ or general relativis- 
tic simulations with simplified microphysics pnjJ22)JSDJ2JJ33 Recently, gravitational 
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waveforms emitted in the rotating core collapse were derived by multidimensional 
simulations in general relativistic frameworks^ 1 '® adopting a finite-temperature nu- 
clear EO£p3 and the electron capture. In their studies, however, the electron capture 
rate was not calculated in a self-consistent manner. Instead, they adopted a simpli- 
fied prescription proposed in Ref. 36) which is based on the result of a spherically 
symmetric simulation. However, it is not clear whether this treatment is justified 
for non-spherical collapse or not. Moreover, they did not take account of emission 
processes of neutrinos. More sophisticated simulations including microphysics are 
required to make accurate predictions of gravitational waveforms. 

The gravitational collapse of massive star is also the primary mechanism of black 
hole formation. Understanding the process of black hole formation is one of the most 
important issues in the theory of the stellar core collapse. A wide variety of recent 
observations have shown that black holes actually exist in the universe (e.g., see 
Ref. l37|) ). and so far, about 20 stellar-mass black holes for which the mass is de- 
termined within a fairly small error have been observed in binary systems of our 
Galaxy and the Large Magellanic Clouds .® The formation of a black hole through 
the gravitational collapse is a highly nonlinear and dynamical phenomenon, and 
therefore, numerical simulation in full general relativity is the unique approach to 
this problem. In spherical symmetry, fully general relativistic simulations of stellar 
core collapse to a black hole have been performed in a state-of-the-art manner, i.e., 
employing realistic EOSs, implementing microphysical processes, and the Boltzmann 
transfer of neutrinos ®>® j n ^ e multidimensional case, by contrast, simulations 
only with simplified microphysics have been performed.®'®'® Because multi- 
dimensional effects such as rotation and convection are likely to play an important 
role, multidimensional simulations in full general relativity employing a realistic EOS 
and detailed microphysics are necessary for clarifying the formation process of black 
holes. 

Furthermore, recent observations®'®'® have found the spectroscopic connec- 
tions between several SNe and long gamma-ray bursts (GRBs) and clarified that 
some of long GRBs are associated with the collapse of massive stars. Supported by 
these observations, the collapsar model 4 ^' is currently one of promising models for 
the central engine of GRBs. In this model, one assumes that the central engine of 
the long GRBs is composed of a rotating black hole and a hot, massive accretion 
disk. Such a system may be formed as a result of the collapse of rapidly rotating 
massive core. In this model, one of the promising processes of the energy-deposition 
to form a GRB fireball is the pair annihilation of neutrinos emitted from the hot, 
massive disk (y e + v e — > e~ + e + ). The collapsar model requires the progenitor core 
to be rotating rapidly enough that the massive accretion disk can be formed around 
the black hole. Recent general relativistic numerical analyses have shown that if a 
progenitor of the collapse is massive and the angular momentum is large enough, 
a black hole surrounded by a massive disk will be formed.®'®'® However, the 
formation mechanism of such system has not been clarified in detail. These also 
enhance the importance of exploring the stellar core collapse to a black hole taking 
account of microphysical processes. 

As reviewed above, multidimensional simulations of stellar collapse in full gen- 
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eral relativity including microphysics is currently one of the most required subjects 
in theoretical astrophysics. However, there has been no multidimensional code in 
full general relativity that self-consistently includes microphysics such as realistic 
EOS, electron capture, and neutrino emission. There have only existed fully general 
relativistic codes in spherical symmetry^ '^S 1 ''^' or Newtonian codes in multidimen- 
sion.^'^'^ We have developed a fully general relativistic multidimensional code 
including a finite-temperature nuclear EOS, self-consistent treatment of the electron 
capture, and a simplified treatment of neutrino emission for the first time. In this 
code, by contrast with the previous onesp^'^ the electron capture rate is treated in 
a self-consistent manner and the neutrino cooling is taken into account for the first 
time. Because it is not currently feasible to fully solve the neutrino transfer equa- 
tions in the framework of general relativity in multidimension because of restrictions 
of computational resources, it will be reasonable to take some approximation for the 
transfer equations at the current status. In this paper, the so-called neutrino leakage 
scheme is adopted as an approximate treatment of neutrino cooling, and a general 
relativistic version of the leakage scheme is developed. 

1.2. The leakage schemes 

The neutrino leakage schemes^ Mi Mi MS Mi Mi Mi Mi as an approximate method 
for the neutrino cooling have a well-established history (e.g. Ref. |57|) ). The basic 
concept of the original neutrino leakage schemes^ 1 is to treat the following two 
regions in the system separately: one is the region where the diffusion timescale of 
neutrinos is longer than the dynamical timescale, and hence, neutrinos are 'trapped' 
(neutrino-trapped region); the other is the region where the diffusion timescale is 
shorter than the dynamical timescale, and hence, neutrinos stream out freely out 
of the system (free-streaming region). The idea of treating the diffusion region sep- 
arately has been applied to more advanced methods for the neutrino transfer (see 
e.g., Ref. |59|) and references therein). 

Then, electron neutrinos and electron anti-neutrinos in the neutrino-trapped 
region are assumed to be in the /3-equilibrium state. The net local rates of lepton- 
number and energy exchange with matters are set to be zero in the neutrino-trapped 
region. To treat diffusive emission of neutrinos leaking out of the neutrino-trapped 
region, simple phenomenological source terms based on the diffusion theory are in- 
troduced.^'^ In the free-streaming region, on the other hand, it is assumed that 
neutrinos escape from the system without interacting with matters. Therefore, neu- 
trinos carry the lepton number and the energy according to the local weak-interaction 
rates. Note that the neutrino fractions are not solved in the original version of the 
leakage scheme: Only the total lepton fraction (from which the neutrino fractions are 
calculated under the /3-equilibrium condition) is necessary in the neutrino-trapped 
region, and the neutrino fractions are set to be zero in the free-streaming region. 
As a result, neutrino quantities and the electron fraction are discontinuous at the 
boundary the neutrino-trapped and free-streaming regions. 

The boundary was given by hand as a single 'neutrino-trapping' density (ptrap) 
without calculating the optical depths of neutrinos in the previous studies.'SD :E3Ji ,[MJ ,[55]i ,[58]i 
However, the location at which the neutrino trapping occurs in fact depends strongly 
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on the neutrino energies (E v ) &2r^ ptrap oc E u ~ 3 , and hence, there are different 
neutrino-trapping densities for different neutrino energies. In the previous leakage 
schemes JSDJS3JS3>EDJSHJ Q n the other hand, all neutrinos were emitted in one moment 
irrespective of their energy. Consequently in the case of the so-called neutrino burst 
emission (e.g., Ref. I60p ). for example, the duration in which the neutrinos are emit- 
ted was shortened and the peak luminosity at the burst was overestimatedP^ 1 '^ 1 '^ 1 
The dependence of the neutrino-trapping densities and the neutrino diffusion rates on 
the neutrino energies are approximately taken into account in the recent simulations 
of mergers of binary neutron starP^'^* However, the lepton-number conservation 
equations for neutrinos are not solved,® which is important to estimate the phase 
space blocking due to neutrinos. 

Transfer equations of neutrinos are not solved in the leakage schemes. There- 
fore, the leakage schemes cannot treat non-local interactions among the neutrinos 
and matters. For example, the so-called neutrino heating^ 1 and the neutrino pair 
annihilation cannot be treated in the leakage scheme. Nevertheless, we believe a 
detailed general relativistic leakage scheme presented in this paper to be a valuable 
approach because even by this approximated approach it is possible to incorporate 
the effects of neutrinos semi-quantitatively as shown in this paper. Also, the neu- 
trino leakage scheme is an appropriate method for studying a number of phenomena 
for which the neutrino heating and neutrino transfer are expected to be not very im- 
portant, e.g., prompt formation of a black hole and compact binary mergers. Both 
of these phenomena are the targets of the present code. 

A first attempt towards a general relativistic leakage scheme was done in the 
previous study. 61 ' In that study, not the region of the system but the energy mo- 
mentum tensor of neutrinos was decomposed into two parts; 'trapped-neutrino' and 
'streaming-neutrino' parts. However the source terms of hydrodynamic and lepton- 
number-conservation equations were determined using the single neutrino-trapping 
density as in the case of the previous leakage schemes. In this paper, we develop a 
new code implementing the microphysical processes in the general relativistic frame- 
work based on the previous study. 61 ' As an application of the code, we perform 
simulations of stellar core collapse. 

A lot of improved ingredients are installed into the present code: (1) The de- 
pendence of the neutrino diffusion rates on the neutrino energies are approximately 
taken into account following the recent study^ with detailed cross sections, instead 
of adopting the single neutrino-trapping density (see Appendix C). (2) The lepton- 
number conservation equations for neutrinos are solved to calculate self-consistently 
the chemical potentials of neutrinos. Then, the blocking effects due to the pres- 
ence of neutrinos and the /3-equilibrium condition can be taken into account more 
accurately (see $3]). (3) A stable explicit method for solving the equations of hy- 
drodynamics, the lepton-number conservations, and neutrinos are developed. Such 
a special procedure is necessary because the characteristic timescale of the weak- 
interaction processes (hereafter referred to as the WP timescale t wp ~ |Y^/YU) is 
much shorter than the dynamical timescale t j vn in hot, dense matter regions 
Note that in the previous leakage schemes EDMKI33>[IIKI5H1> the /3-equilibrium was as- 
sumed to be achieved in such regions (i.e. Y e = 0) and no such special treatments 
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are required. See ^2] for further discussions and jj^for details of the method. (4) The 
electron capture rate are calculated in a detailed manner®} including effects of the 
so-called thermal unblocking^ (see Appendix A). 

The paper is organized as follows. First, issues in implementation of weak inter- 
actions and neutrino cooling in full general relativistic simulation is briefly summa- 
rized in £}2j Then, framework of the general relativistic leakage scheme is described 
in detail in § [3j In § HJ EOSs employed in this paper are described in some details. 
Basic equations and numerical methods of the simulations are described in § [5j Nu- 
merical results obtained in this paper are shown in § [6j We devote § [7] to a summary 
and discussions. In appendices, details of the microphysics adopted in the present 
paper are summarized for the purpose of convenience. Throughout the paper, the 
geometrical unit c = G = 1 is used otherwise stated. 

§2. Issues in implementation of weak interactions and neutrino cooling 
in fully general relativistic simulation 

Because the characteristic timescale of the weak-interaction processes (the WP 
timescale t wp ~ |Y e /l^|) is much shorter than the dynamical timescale tdyn m hot 
dense mattersj^'^ the explicit numerical treatments of the weak interactions are 
computationally expensive in simple methods, as noted in the previous pioneering 
work by BruennP*3 A very short timestep (At < t wp <C idyn) will be required to 
solve the equations explicitly. 

The net rates of lepton-number and energy exchanges between matters and neu- 
trinos may not be large, and consequently, an effective timescale may not be as short 
as the dynamical timescale. However, this does not immediately imply that one can 
solve the equations explicitly without employing any prescription. For example, the 
achievement of ^-equilibrium, where Y e = is the consequence of the cancellation 
of two very large weak interaction processes (the electron and the electron-neutrino 
captures, see Eq. (|3-20p ) and of the action of the phase space blocking. Note that 
the weak interaction processes depend enormously both on the temperature and the 
lepton chemical potentials. Therefore, small error in the evaluation of the tempera- 
ture and a small deviation from the /3-equilibrium due to small error in calculation 
of the lepton chemical potentials will result in huge error. Then, stiff source terms 
appear and explicit numerical evolution often becomes unstable. Indeed, we found 
that a straightforward, explicit solution of the equations did not work. 

In the following of this section, we describe issues of implementation of weak 
interactions and neutrino cooling into the hydrodynamic equations in the conserva- 
tive schemes in fully general relativistic simulations. Fiest, we illustrate that in the 
Newtonian framework, the equations may be solved implicitly in a relatively sim- 
ple manne ^i mm mm mm mm ( see also Refs. EH) and and references 
therein). The equations of hydrodynamics, lepton-number conservations, and neu- 
trino processes are schematically written as, 



p = 0, 

Vi = S Vi (p,Y e ,T,Q u ), 



(2-1) 
(2-2) 
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Y e = S Ye (P,Y e ,T,Q u ), (2-3) 
e = S e (p,Y e ,T,Q u ), (2-4) 
Qv = Sq u (p, Y e , T, Q u ), (2-5) 

where p is the rest-mass density, Vi is the velocity, Y e is the electron fraction, e is 
the (internal) energy of matter, T is the temperature, and Q u stands for the relevant 
neutrino quantities. We here omit the transport terms. S's in the right-hand side 
stand for the relevant source terms. Comparing the quantities in the left-hand-side 
and the argument quantities in the source terms, only the relation between e and T 
is nontrivial. Usually, EOSs employed in the simulation are tabularized, and one di- 
mensional search over the EOS table is required to solve them. Due to the relatively 
simple relations between the quantities to be evolved and the argument quantities, 
the above equations may be solved implicitly in a straightforward (although compli- 
cated) manner. 

In the relativistic framework, the situation becomes much more complicated 
in conservative schemes, because the Lorentz factor (i -1 ) is coupled with rest-mass 
density and the energy density (see Eqs. (|5-12p and (|5-16p where w = au 1 " is used 
instead of r), and because the specific enthalpy (h = h(p, Y e ,T)) is coupled with the 
momentum (see Eq. (|5-14|) ). 

It should be addressed that the previous fully general relativistic works in the 
spherical symmetry^ are based on the so-called Misner-Sharp coordinates.^ 1 
There are no such complicated couplings in these Lagrangian coordinates. Accord- 
ingly, the equations may be solved essentially in the same manner as in the Newtonian 
framework. Because no such simple Lagrangian coordinates are known in the mul- 
tidimensional case, the complicated couplings inevitably appear in the relativistic 
framework. 

Omitting the factors associated with the geometric variables (which are usually 
updated before solving the hydrodynamics equations) and the transport terms, the 
equations to be solved in the general relativistic framework are schematically written 

as, 

p*(p,r) = 0, (2-6) 
Ui(ui,h) = Ui(ui,p,Y e ,T) = Su t (p,Y e ,T,Q u ,r), (2-7) 

Y e = S Ye (p,Ye,T,Q u ,r), (2-8) 
k(p, Y e ,T, r) = Sg(p, Y e , T, Q u , r), (2-9) 
Q u = S Qu (p,Y e ,T,Q u ,r), (2-10) 

where p* is a weighted density, u a is a weighted four velocity, e is a weighted energy 
density (see § 15.21 for the definition of these variables). The Lorentz factor is calcu- 
lated by solving the normalization condition u a u a = — 1, which is rather complicated 
nonlinear equation schematically written as 

/normalization {U{ , — /normalization (^i j Pi Yet T, ij — 0. (2'11) 

The accurate calculation of the Lorentz factor and the accurate solution of the nor- 
malization condition are very important in the numerical relativistic hydrodynamics. 



8 



Y. Sekiguchi 



Now, it is obvious that the argument quantities in the source terms are not 
simply related with the evolved quantities in the left-hand-side of Eqs. (|2-6[) — ()2- lip . 
Solving the equations implicitly is not as straightforward as in the Newtonian case 
and no successful formulations have been developed. Moreover it might be not clear 
whether a convergent solution can be stably obtained numerically or not, because 
they are simultaneous nonlinear equations. Therefore, it may be not a poor choice 
to adopt an alternative approach in which the equations are solved explicitly with 
some approximations as described in the next sectiorPl. 

§3. General relativistic neutrino leakage scheme 

In this section, we describe a method for approximately solving hydrodynamic 
equations coupled with neutrino radiation in an explicit manner. As described in the 
previous section, because of the relation i wp <C tdyn in the hot dense matter regions, 
the source terms in the equations become too stiff for the equations to be solved 
explicitly in the straightforward manner. The characteristic timescale of leakage of 
neutrinos from the system ii ea k) by contrast, is much longer than t wp in the hot dense 
matter region. Rather, ii ea k ~ L/c ~ tdyn where L is the characteristic length scale 
of the system. On the other hand, ii ea k is comparable to t wp in the free-streaming 
regions but t wp is longer than or comparable with idyn there. All these facts imply 
that the WP timescale does not directly determine the evolution of the system but the 
leakage timescale does. Using this fact, we approximate some of original equations 
and reformulate them so that the source terms are to be characterized by the leakage 
timescale ti ea k- 

3.1. Decomposition of neutrino energy-momentum tensor 

The basic equations of general relativistic hydrodynamics with neutrinos are 

Va(T Total ) a = Va [ (T F)a + ^a] = ^ (g.-Q 

where (T Total ) Q , ( g is the total energy-momentum tensor, and (T F ) a p and {T L ') a p are 
the energy- momentum tensor of fluids and neutrinos, respectively. Equation (|3-ip 
can be written as 

V a {T ¥ )1 = Qp, (3-2) 
V Q (T^ = -Qp, (3-3) 

where the source term Q a is regarded as the local production rate of neutrinos 
through the weak processes. 

Now, the problem is that the source term Q a becomes too stiff to solve explicitly 
in hot dense matter regions where i wp <C idyn- To overcome this situation, the 
following procedures are adopted. 

First, it is assumed that the energy- momentum tensor of neutrinos are be decom- 
posed into 'trapped-neutrino' ((T" / ' T ) Q!/ g) and 'streaming-neutrino' ((T u ' S ) a p) parts 

*^ It should be stated that the implicit schemes are also approximated ones because a short WP 
timescale associated with the weak interaction is not fully resolved. 
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as,® 

{T V U = ( TU ' T U + {T^U- (3-4) 
Here, the trapped-neutrinos phenomenologically represent neutrinos which interact 
sufficiently frequently with matter and are thermalized. On the other hand, the 
streaming-neutrino part describes a phenomenological flow of neutrinos streaming 
out of the system^ (see also Ref. I76|) where a more sophisticate method in terms of 
the distribution function is adopted in the Newtonian framework). 

Second, the locally produced neutrinos are assumed to leak out to be the streaming- 
neutrinos with a leakage rate Qj^ ak : 

V?{T v > B t = Q 1 ^- (3-5) 

Then, the equation of the trapped-neutrino part is 

Vp(T^)i = Q a -Q 1 ^. (3-6) 

Third, the trapped-neutrino part is combined with the fluid part as 

T a p = (T F ) a/3 + (T»> T ) a p, (3-7) 

and Eqs. (|3-2[) and (|3-6|) are combined to give 

V/^ = -Q 1 ^. (3-8) 



Thus, the equations to be solved are changed from Eqs. (|3-2p and f|3-3|) to Eqs. f|3-8j) 
and (|3-5p . Note that the new equations only include the source term Q^ 3 ^ which is 
characterized by the leakage timescale ii ea k- Definition of Q^ ak will be given in § 13.31 
The energy- momentum tensor of the fluid and trapped-neutrino parts (T a p) is 
treated as that of the perfect fluid, 

T a p = + pe + P)u a up + Pg aP , (3-9) 

where p and u a are the rest mass density and the 4-velocity. The specific inter- 
nal energy density (e) and the pressure (P) are the sum of contributions from the 
baryons (free protons, free neutrons, a-particles, and heavy nuclei), leptons (elec- 
trons, positrons, and trapped-neutrinos), and the photons as, 

P = P B +P e + P u + P ph , (3-10) 
e = £ B + e e + e v +e ph , (3-11) 

where subscripts 7 P>\ 'e\ 'ph', and V denote the components of the baryons, elec- 
trons and positrons, photons, and trapped-neutrinos, respectively. 

The streaming-neutrino part, on the other hand, is set to be a general form of 

(T u ' S ) a p = En a np + F a n p + F p n a + P a(5 , (3-12) 

where F a n a = P a /3n a = 0. In order to close the system, we need an explicit expres- 
sion of P a p. In this paper, we adopt a simple form P a p = xEla/3 with x = 1/3- 
This approximation may work well in high density regions but will violate in low 
density regions. However, the violation will not affect the dynamics because the 
total amount of streaming-neutrinos emitted in low density regions will be small. Of 
course, a more sophisticated treatment will be necessary in a future study. 
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3.2. The lepton-number conservation equations 

The conservation equations of the lepton fractions are written schematically as 

(3-13) 

(3-14) 

(3-15) 

(3-16) 



dYe _ 




dt 


= -7e, 


dY ue 




dt 


= lue, 


dY Pe 




dt 


= lue, 


dY 

<* 1 ux 




dt 


= lux 



lue 


local 


_ ~ leak 


lue 


lue ■ 


lue 


— - local 


_ ~ leak 


lue 


lue • 


lux 


— ^local 


leak 


1 ux 


lux 



where Y e , Y ue , Y ue , and Y ux denote the electron fraction, the electron neutrino frac- 
tion, the electron anti-neutrino fraction, and /j, and r neutrino and anti-neutrino 
fractions, respectively. We note that in the previous simulations based on the leak- 
age schemes ED JSZ) ,EH]),|63} neiUT i no fractions were not solved. 

The source terms of neutrino fractions can be written, on the basis of the present 
leakage scheme, as 

(3-17) 
(3-18) 
(3-19) 

where <-y local 'g and ^ leak 's are the local production and the leakage rates of each neu- 
trino, respectively (see § I3.3j> . Note that only the trapped- neutrinos are responsible 
for the neutrino fractions. Assuming that the trapped neutrinos are thermalized and 
the distribution function is the equilibrium Fermi-Dirac one, the chemical potentials 
of neutrinos can be calculated from the neutrino fractions. Then the thermodynam- 
ical quantities of neutrinos can be also calculated. 

The source term of the electron fraction conservation is 

7e = 7 1 °e Cal -7pe Cal - (3-20) 

Because ^l ocal 'g are characterized by the WP timescale i wp , some procedures are 
necessary to solve the lepton conservation equations explicitly. The following simple 
procedures are proposed to solve the equations stably. 

First, in each timestep n, the conservation equation of the total lepton fraction 
(Y, = Y e -Y m + Y K ), 

f " -* < 3 ' 21 > 

is solved together with the conservation equation of Y ux , Eq. (|3-16p . in advance of 
solving whole of the lepton conservation equations (Eqs. (|3-13p - (|3-16p ). Note that 
the source term 'ji — l]^^ — Tp^ a ^ ^ 

characterized by the leakage timescale ti ea k so 
that this equation can be solved explicitly in the hydrodynamic timescale. Then, 
assuming that the /3-equilibrium is achieved, values of the lepton fractions in the 
/3-equilibrium (Y e , Y„ e , and Yp e ) are calculated from the evolved Y/. 
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Second, regarding Yue and Yp e as the maximum allowed values of the neutrino 
fractions in the next timestep n + 1, the source terms are limited so that Y^'s in 
the timestep n + 1 cannot exceed Yy 's. Then, the whole of the lepton conservation 
equations (Eqs. (|3-13p - (|3-16p ) are solved explicitly using these limiters. 

Third, the following conditions are checked 

flp + fl e < fi n + /i„ e , (3-22) 

Hn ~ He < Hp + Hue, (3-23) 

where Hp, Hn, He, Hue and Hue are the chemical potentials of protons, neutrons, 
electrons, electron neutrinos, and electron anti-neutrinos, respectively. If both con- 
ditions are satisfied, the values of the lepton fractions in the timestep n + 1 are set to 
be those in the /3-equilibrium value; Ye , Y{? e , and Yp e . On the other hand, if either 
or both conditions are not satisfied, the lepton fractions in the timestep n + 1 is set 
to be those obtained by solving whole of the lepton-number conservation equations. 

A limiter for the evolution of Y ux may be also necessary in the case where the pair 
processes are dominant, for example, in the simulations for collapse of population III 
stellar core. In this case, the value of Y ux at the pair equilibrium (i.e. at Hvx = 0), 
Yyx U may be used to limit the source term. 

3.3. Definition of leakage rates 

In this subsection the definitions of the leakage rates Q^ ak and <yjf ak are pre- 
sented. Because Q|f ak may be regarded as the emissivity of neutrinos measured in 
the fluid rest frame, Q^ ak is defined a#B 

gteak = gtaak^ (3 . 24) 



The leakage rates Q^ 3 ^ and ^f ak are assumed to satisfy the following properties. 

1. The leakage rates approach the local rates Q^° cal and 

7 locaI in 

the low density, 

transparent region. 

2. The leakage rates approach the diffusion rates Q„ and 7^ iff in the high density, 
opaque region. 

3. The above two limits should be connected smoothly. 

Here, the local rates can be calculated based on the theory of weak interactions and 
the diffusion rates can be determined based on the diffusion theory (see appendices 
for the definition of the local and diffusion rates adopted in this paper). There will be 
several prescriptions to satisfy the requirement (ii^pU'ES j n paper, the leakage 
rates are defined by 

g|f* = (1 - e - bT »)Qf s + e- hT »Q x ° c& \ (3-25) 

7 leak = ^ _ e -&r„) 7 diff + g-br^loca^ (3.^ 

where t v is the optical depth of neutrinos and b is a parameter which is typically 
set as b~ x = 2/3. The optical depth can be computed from the cross sections in a 
standard manner ES^EU 



12 



Y. Sekiguchi 



In the present implementation, it is not necessary to artificially divide the system 
into neutrino-trapped and free-streaming regions by the single neutrino-trapping 
density. Therefore there is no discontinuous boundary which existed in the previous 
leakage schemes jSDJSDJUJ 

As the local production reactions of neutrinos, the electron and positron cap- 
tures® (7^ and 7^); the electron-positron pair annihilation^ (7*^ for electron- 
type neutrinos and for the other type), the plasmon decays 63 ' (7^^ and 7^.^,), 
and the Bremsstrahlung processes^} (7^ r p™ s and 7^. r p™ s ) are considered in this pa- 
per. Then, the local rates for the neutrino fractions are 

^local ec , pair plas Brems (3-27) 
ive lue ' iu e u e ' lu e u e ' lu e i/ e i \" ) 

local pc pair plas Brems /g^gN 

7 i° cal = 7 p^ + 7^ + 7^ r p ms - (3-29) 
Similarly, the local neutrino energy emission rate Q}? cal is given by 

local /Qcc , /^)PC , o /^pair . ^plas . ^BremsN 



I'cl'c 



+ 4 {Ql% + Ql% + QfX 5 ) ■ (3-30) 

The explicit forms of the local rates in Eqs. (|3-27p - (|3-30p will be found in Appendices 
A and B. 

We follow the recent work by Rosswog and Liebenddrfer 62 ' for the diffusive 
neutrino emission rates 7 ^ ifT and Qf S in Eqs (EE25|I and flESD- The explicit forms 
of 7^ and Q„ lS are presented in Appendix C. 



§4. Equation of state 



In this section we summarize details of EOSs adopted in our current code. 
4.1. Baryons 

In the present version of our code, we employ an EOS by Shen et aL,™ which 
is derived by the relativistic mean field theory®} based on the relativistic Briickner- 
Hartree-Fock theory!® The so-called parameter set TMl®) is adopted to reproduce 
characteristic properties of heavy nuclei. The maximum mass of a cold spherical neu- 
tron star in this EOS is much larger than the canonical neutron star mass ~ 1.4M Q 
as ~ 2.2M & ?^ The framework of the relativistic mean field theory is extended 
with the Thomas-Fermi spherical cell model approximation to describe not only the 
homogeneous matter but also an inhomogeneous one. 

The thermodynamical quantities of dense matter at various sets of (p, Y p , T) are 
calculated to construct the numerical data table for simulation. The table covers a 
wide range of density 10 51 -10 15 ' 4 g/cm 3 , electron fraction 0.0-0.56, and temperature 
0-100 MeV, which are required for supernova simulation. It should be noted that the 
causality is guaranteed to be satisfied in this framework, whereas the sound velocity 
sometimes exceeds the velocity of the light in the non-relativistic framework, e.g., 
in the EOS by Lattimer and SwestyPZl This is one of the benefits of the relativistic 
EOS. 
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Although we employ the nuclear EOS by Shen et al. in this work, it is easy to 
replace the EOS. In the future we plan to implement other EOSs such as a hyperonic 
matter EOSP 

Because the table of the original EOS by Shen et al. does not include the 
thermodynamical quantities of the leptons (electrons, positrons, and neutrinos if 
necessary) and photons, one has to consistently include them to the table. 

4.2. Electrons and Positrons 

To consistently calculate the pressure and the internal energy of the electron and 
positron, the charge neutrality condition Y p = Y e should be solved to determine the 
electron chemical potential \x e for each value of the baryon rest-mass density p and 
the temperature T in the EOS table. Namely, it is required to solve the equation 

n e (p e ,T) = n_ - n+ = — (4-1) 

m u 

in terms of p e for given values of p, T, and Y e (= Y p ). Here, m u = 931.49432 MeV is 
the atomic mass unit, and n_ and n+ are the total number densities (i.e., including 
the electron-positron pair) of the electrons and positrons, respectively. 

Assuming that the electrons obey the Fermi-Dirac distribution (which is derived 
under the assumption of thermodynamic equilibrium), the number density (n_), the 
pressure (-P-), and the internal energy density (ti-) of the electron are written as® 

n = j_ r (4 .2) 

ttWo exp[-77 e + e~A B T] + l' 1 > 

J_ f°° P Hdi/dp)dp 
- tt 2 Wo exphTfe + e/fcjjTl + l' 1 ' 

1 [°° P 2 edp 
U - ttWo ew [- Ve + ~ e /k B T] + l 1 ' 

Here h, fee, and rj e = p e /ksT are the Planck's constant, the Boltzmann's constant 
and the so-called degeneracy parameter. e(p) = yj m\c^ + p 2 — m e c 2 is the kinetic 
energy of a free electron. If we choose the zero point of our energy scale for electrons 
at e = 0, we have to assign a total energy of e + = sj m 2 c 4: + p 2 + m e c 2 to a free 
positron.® Thus the number density (n+), the pressure (P+), and the internal 
energy density (u + ) of positrons are given by® 

(4-5) 
(4-6) 
(4-7) 



n+ = 




poo 


p 2 dp 




ir 2 h? . 


Jo ex P 


[-V+ + e+/k B T\ 


+ 1 


P+ = 




POO 


p 3 (de + /dp)dp 




ir 2 h? . 


Jo ex P 


[-V+ + h/kEsT] 


+ 1 


u + = 




POO 


p 2 (e + 2m e c 2 )dp 




n 2 h 3 . 


Jo ex P 


[-ri + + e+/k B T] 


+ 1 



where r/ + = 



—7] e is the degeneracy parameter of the positrons. 
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4.3. Photons 

The pressure and the specific internal energy density of photons are given by 



Pr = 



a r T A 



a r T A 



e r = 



(4-8) 



3 p 

where a r is the radiation constant a r = (ir 2 k A B ) / (15c 3 h 3 ) and c is the velocity of light. 

4.4. Trapped neutrinos 

In this paper, the trapped-neutrinos are assumed to interact sufficiently fre- 
quently with matter that be thermalized. Therefore they are described as ideal 
Fermi gases with the matter temperature. Then, from the neutrino fractions Y u , the 
chemical potentials of neutrinos are calculated by solving 



Y u = Y u (p u ,T). 



(4-9) 



Using the chemical potentials, p u , and the matter temperature, the pressure and the 
internal energy of the trapped-neutrinos are calculated in the same manner as for 
electrons. 

4.5. The sound velocity 

In the high-resolution shock-capturing scheme for hydrodynamics, we in general 
need to evaluate the sound velocity c s , 



1 



dP 



dp 



+ 



P dP 

p de 



(4-10) 



The derivatives of the pressure are calculated by 



dP_ 

dp 

dP_ 

de 



- E 



=B,e,r,i/ 



m 

dp 



m 

dT 



E 



\j=B,e,r,v 



dp 






(4-11) 
(4-12) 



where '-B', 'e', 'p/i' and V in the sum denote the baryon, electron, photons, and 
neutrino quantities. 

The derivatives for the baryon parts are evaluated by taking a finite difference 
of the table data. On the other hand, the derivatives for the electron parts can be 
evaluated semi-analytically. Because there are in general the phase transition regions 
in an EOS table for baryons and moreover the EOS may contain some non-smooth 
spiky structures, careful treatments are necessary when evaluating the derivatives of 
thermodynamical quantities. In the present EOS table, the derivatives are carefully 
evaluated so that there are no spiky behaviors in the resulting sound velocities. 
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§5. Basic equations and Numerical methods 

5.1. Einstein's equation and gauge conditions 

The standard variables in the 3+1 decomposition are the three-dimensional met- 
ric 7ij and the extrinsic curvature K{j on the three-dimensional hypersurface® de- 
fined by 

Ip-v = 9^u + n^n u , (5-1) 

KfW = — Ij^-nJ 'fj,i> > 

where g^ v is the spacetime metric, is the unit normal to a three-dimensional 
hypersurface, and £ n is the Lie derivative with respect to the unit normal n^. Then 
we can write the line element in the form 

ds 2 = -a 2 dt 2 + jijidx 1 + /3 l dt)(dx j + ftdt), (5-3) 

where a and j3 l are the lapse function and the shift vector which describe the gauge 
degree of freedom. 

In the BSSN reformulation ESJJHD the spatial metric jij is conformally decom- 
posed as jij = e^jij where the condition det(7y) = 1 is imposed for the conformal 
metric jij. From this condition, the conformal factor is written as <p = i In 7 and 
7 = det(7y). The extrinsic curvature is decomposed into the trace part K 
and the traceless part Aij as Kij = Aij + (1/3) jij K . The traceless part Aij is 
conformally decomposed as = e^Aij. Thus the fundamental quantities for the 
evolution equation are now split into (j>,jij, K, and Ay. Furthermore, the auxiliary 
variable Fi = S^ k d k jij is introduced in the BSSN reformulation.® 

The basic equations to be solved are 

d t - /3 k d k ) <f> = i {-aK + d k f3 k ) , (5-4) 



P k d k ) j i3 = -2aAij + j ik dj f3 k + j ]k dip k - ^j i3 d k (3 k , (5-5) 



dt 

dt ~ P k d k ) K = -D k D k a + a 



d t - /3 k d k ) = 



Aijtfi + -if 2 



3 

1 

3 



+ 4vra ( Ph + S) , (5-6) 



a [Rij - \e A(t> jijR \ - {c>iD 3 a - ^e 4 %jD k D k a) 
+a (KAij - 2A ik A k j) + A ik djP k + A jk d^ k - \A l3 d k p k 



3 

-8vra ( e'^Sij - U^s) , (5-7) 



(d t - P k d k ^ Fi = -lforaj, 

+2a lf kj djA ik + A ik djf k i - X -Ai [ dih 3l + 61*^0 - 
+5 ik [-2A\jd k a+ (dk^dihij 
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+d k (ludrf 1 + jjAP 1 - | , (5-8) 

where ^R, ^R ij: and A are the Ricci scalar, the Ricci tensor, and the covariant 
derivative associated with three-dimensional metric 7™, respectively. The matter 
source terms, p h = (T Total )^n a n^, £ = -(T Total )^ 7iQ n^, and 5 - = (T Total ) a ^ ialj p, 
are the projections of the stress-energy tensor with respect to and 7^, and 
S = -, r S".'. 

We assume the axial symmetry of the spacetime and the so-called Cartoon 
method®}'^ is adopted to avoid problems around the coordinate singularities of 
the cylindrical coordinates. Except for this, the numerical schemes for solving the 
Einstein's equation are essentially the same as those in Ref. [90]). We use 4th-order fi- 
nite difference scheme in the spatial direction and the 3rd-order Runge-Kutta scheme 
in the time integration. The advection terms such as f3 l di(fi are evaluated by a 4th- 
order upwind scheme. 

As the gauge conditions for the lapse, we use the so-called 1 + log slicing® 1 

(dt - £p)a = -2Ka. (5-9) 

It is known that the 1 + log slicing enables to perform a long term evolution of 
neutron stars as well as has strong singularity avoidance properties in the black hole 
spacetime. 

The shift vector is determined by solving the following dynamical equation^ 

d t p k = f\F l + Atd t F l ). (5-10) 
Here the second term in the right-hand side is necessary for numerical stability.^ 1 

5.2. The hydrodynamic equations in leakage scheme 

The basic equations for general relativistic hydrodynamics in our leakage scheme 
are the continuity equation, the lepton-number conservation equations, and the local 
conservation equation of the energy-momentum. We assume the axial symmetry of 
the spacetime and the hydrodynamics equations are solved in the cylindrical coordi- 
nates (w, cp, z) where zu = \J x 1 + y 2 . In the axisymmetric case, the hydrodynamics 
equations should be written in the cylindrical coordinate. On the other hand, in the 
Cartoon method)® Einstein's equation are solved in the y = plane for which 
x = ru, u m = u x , u v = xu y , and the other similar relations hold for vector and ten- 
sor quantities. Taking into these facts, the hydrodynamic equations may be written 
using the Cartesian coordinates replacing (w, (p) by (x, y). In the following, we write 
down explicit forms of the equations for the purpose of convenience. Numerical tests 
for basic parts of the code of solving the hydrodynamics equations are extensively 
performed in Ref. 89). The equations are solved using the third-order high-resolution 
central scheme of Kurganov and TadmorP^'EU 

5.2.1. The Continuity and lepton-number conservation equations 
The continuity equation for the baryon rest mass is 



V a (pu a ) = 0. 



(5-11) 
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As fundamental variables for numerical simulations, the following quantities are in- 
troduced: = pwe 6 ^ and v l = u l ju l where w = au l . Then, the continuity equation 
is written as 



dt(p*) + -d x (p*v x ) + d z (p*v 2 



0. 



(5-12) 



Using the continuity equation, the lepton-number conservation equations (|3T3p 
()3-16p are written as 



d t (p*Y L ) + ~d x {p*Y L v x ) + d z (p*Y L v z ) = p^ L , 
x 



(5-13) 



where Yl and are abbreviated expressions of the lepton fractions and the source 
terms. 

5.2.2. Energy-momentum conservation 

As fundamental variables for numerical simulations, we define the quantities 
Ui = hui and e = hw — P(pw)~ 1 . Then, the Euler equation jfVpT a = — 7f Qa ak , 
and the energy equation n^V^T? = — n Q Qj^ ak can be written as 



d t (p*u A ) + -d x 
x 



-p* 



x ^p*uav x + 



+ d z 



P*UAV" 



+ Pae^S 2 



whdAa — UidAfi 1 + 



ae 



2wh 



-u k uid A l '" 



2ah(w 2 - 1) 



d A (f) 



+Pd A (ae 6(t> ) + yr *~ y " ' - — >" A - ae^QT k , (5-14) 



(p*u y yy + Pae^)6 A _ 

(Jits \a£ A 

X 

d t (p*u y ) + ^d x {x 2 p,u y v y ) + d z {p,u y v z ) = -ae^Q 1 ^, 
d t {p*e) + ^8 X [x [p*v x e + Pe &<t> {v x + /3 X )} 



+ d z 



p*v e + 



v z + P z ) 



ae^PK + -^-u k uiK kl - p^W-jO - ae b<p Q l ™ K n a , 



leat: a 



u l h 

where the subscript A denotes x or z component. 

The evolution equation of streaming-neutrinos V 0(T u ' S )a = Q„ ak gives 



(5-15) 



(5-16) 



dt(E) + -d x x(aF x -p x E) +8 Z (aF z -(3 z E) 



aEK 



F k d k a + ae 6 *Q l ™ k n a , 



(5-17) 



d t {F A ) + -d x 
x 



x ( -aE8 x A - L3 X F A 



+ d z 



-aEb z A - I3 Z F A 



-Ed aol + F k d A P + 2aEd A <l> + 



k , „ , (E/3 - F y f3 y )5 A _ ,„ v .,],,, k 



d t (F y ) ~ ~ 2 d X 



x 2 (3 x F y 



d z 



+ ae b *QT k , (5-18) 
(5-19) 



where E = e^E and F{ = e 6 ^^, and the subscript A again denotes x or z compo- 
nent. The closure relation P a p = E^ a p/?> is also substituted. 
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5.3. Recover of (p, Y e /Y\, T) 

The quantities numerically evolved in the relativistic hydrodynamics are the 
conserved quantities (p*, Hi, e) and the lepton fraction Y e or Y\. The argument 
variables, (p, (Y e or Yj), T), of the EOS table, together with the weight factor 
to = \J\ + "f % iuiUj, should be calculated from the conserved quantities at each time 
slice. Note that the electron (Y e ) or lepton fraction (YJ) is readily given by numerical 
evolution at each time slice whereas p, ui, and T are not. This fact requires us to 
find an efficient method for determining w. 

5.3.1. Non- /3-equilibrium case 

In the case that the /3-equilibrium condition is not satisfied, the argument quan- 
tities (p, Y e , T) can be reconstructed from the conserved quantities in the following 
straightforward manner. 

1. Give a trial value of w, referred to as w. Then, one obtains a trial value of the 
rest mass density from p = p*/{we^). 

2. A trial value of the temperature, T, can be obtained by solving the following 
equation: 

e= | l + e + Z)w-^ = e(p,Y e ,f). (5-20) 
\ P J pw 

Here, one dimensional search over the EOS table is required to obtain T. 

3. The next trial value of w is given by w = \J 1 + e -4 ^^*- 3 UiUjh~ 2 , where the 

specific enthalpy was calculated as h = h(p,Y e ,T) in the step 2. 

4. Repeat the procedures (l)-(3) until a required degree of convergence is achieved. 
Convergent solutions of the temperature and w are obtained typically in 10 
iterations. 

5.3.2. The /3-equilibrium case 

On the other hand, in the case that the /3-equilibrium condition is satisfied, one 
has to reconstruct the argument quantities (p,Y e ,T) from the conserved quantities 
and Yi, under the assumption of the /3-equilibrium. In this case, two-dimensional 
recover (Y;,e) =^ (Y e ,T) would be required for a given value of w. A serious 
problem is that in this case, there may be more than one combination of (Y e , T) 
which gives the same values of Y\ and e. Therefore, we have to adopt a different 
method to recover (p, Y e ,T). Under the assumption of the /3-equilibrium, the electron 
fraction is related to the total lepton fraction: Y e = Y e (p,Yi,T). Using this relation, 
the EOS table can be rewritten in terms of the argument variables of (p, Yj, T). 
Then, the same strategy as in the non- /3-equilibrium case can be adopted. Namely, 

1. Give a trial value w. Then one obtains a trial value of the rest mass density. 

2. A trial value of the temperature can be obtained by solving e = e(p, Yj, T), with 
one dimensional search over the EOS table. 

3. The next trial value of w is given by w = \J 1 + c^^iiiUjh^ 2 . 

4. Repeat the procedures (l)-(3) until a required degree of convergence is achieved. 
The electron fraction is given as Y e = Y e (p,Yi,T) in the (new) EOS table. 
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Model 




<P C < 0.0125 


<<P C < 0.025 


< $ c < 0.05 


< <P C < 0.1 


$c > 0.1 


S15 


Ax 


3.26 


1.60 


0.820 


0.414 


0.217 




5 


0.002 


0.002 


0.002 


0.002 


0.002 




N 


444 


668 


924 


1212 


1532 




L (km) 


2330 


2239 


2188 


2124 


2103 


S15 


Ax 


5.10 


2.90 


1.44 


0.760 


0.396 


(low) 


8 


0.002 


0.00215 


0.0023 


0.00245 


0.0026 




N 


316 


444 


636 


828 


1020 




L (km) 


2244 


2151 


2073 


2043 


2000 



Table I. Summary of the regridding procedure. The values of the minimum grid spacing Axo (in 
units of km), the non-uniform-grid factor 8, and the grid number N for each range of $ c — l — a c 
are listed. 

In the case of a simplified or analytic EOS, the Newton-Raphson method may 
be applied to recover the primitive variables. In the case of a tabulated EOS, by con- 
trast, the Newton-Raphson method may not be a good approach because it requires 
derivatives of thermodynamical quantities which in general cannot be calculated 
precisely from a tabulated EOS by the finite differentiating method. 

5.4. Grid Setting 

In numerical simulations, we adopt a nonuniform grid, in which the grid spacing 
is increased as 

dxj + i = (1 + 5)dxj, dzi + i = (1 + 5)dzi (5-21) 

where dxj = x j+i ~ x j-> dzi = zi + \ — z\, and 5 is a constant. In addition, a regridding 
technique^ is adopted to assign a sufficiently large number of grid points inside 
the collapsing core, saving the CPU time efficiently. The regridding is carried out 
whenever the characteristic radius of the collapsing star decreases by a factor of a 
2-3. At each regridding, the minimum grid spacing is decreased by a factor of ~ 2 
while the geometrical factor 8 is unchanged (see Table [TJ) . 

All the quantities on the new grid are calculated using the fifth-order Lagrange 
interpolation. To avoid discarding the matter in the outer region, we also increase 
the grid number at each regridding. For the regridding, we define a relativistic 
gravitational potential <P C = 1 — a c (# c > 0) where a c is the central value of the 
lapse function. Because <P C is approximately proportional to M/R where M and 
R are characteristic mass and radius of the core, < P~ 1 can be used as a measure of 
the characteristic length scale of the stellar core for the regridding. In Table HI we 
summarize the regridding parameters of each level of the grid. 

§6. Results 

As a test problem, we perform a collapse simulation of spherical presupernova 
core and compare the results with those in the previous studies, to see the validity of 
the present code. Most of the following results are not novel astrophysically, but are 
novel in the sense that stellar core collapse can be followed by a multidimensional 
fully general relativistic simulation taking account of microphysical processes. In 
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10 20 

Time after bounce t b [ms] 



Fig. 1. Evolution of the central density p c (upper panel) and the central value of the lapse function 
« c (lower panel). The solid curves are results for the finer grid resolution and the dotted curves 
are results of the coarser grid resolution. 



§ 16. 2| § 16.31 and § 16. 4| we first review the basic features of the collapse dynamics and 
the shock formation, the stall of shock, and convective activities. Then we compare 
our results with those in the previous studies in § 16.51 

6.1. Initial condition 

In this paper, we adopt a recent presupernova model of massive star by Woosley, 
Heger, and Weaver P^J 15Mq model with solar metallicity (hereafter S15 model). 
We follow the dynamical evolution of the central part which constitutes the Fe core 
and the inner part of the Si-shell. We read in the density, the electron fraction, 
the temperature and the velocity {vi) of the original initial data and derive other 
thermodynamical quantities using the EOS table. 

Note that the procedure of remapping the original initial data into the grid 
adopted in the numerical simulations is coordinate-dependent in general relativity. 
In this paper, we read in the original data as a function of the coordinate radius. In 
this case, the baryon rest-mass of the core is slightly larger than the original one, 
because it is defined by 



M* = J p*dx 3 = J p(we 6,p )d 3 x, 

where we^ > 1. 



Full GR simulation with microphysics 



21 



6.2. Core bounce and shock formation 

Figure [T] displays the time evolution of the central rest-mass density p and the 
central value of the lapse function. This figure shows that the stellar core collapse to 
a neutron star can be divided into three phases; the infall phase, the bounce phase, 
and the quasi-static evolution phase (see Refs. l2Tj) and l2l|) for the case of rotational 
collapse). The general feature of the collapse is as follows. 

The infall phase sets in due to gravitational instability of the iron core trig- 
gered by the sudden softening of the EOS, which is associated primarily with the 
electron capture and partially with the photo-dissociation of the heavy nuclei. The 
collapse in an early phase proceeds almost homologously. However, the collapse in 
the central region is accelerated with time because the electron capture reduces the 
degenerate pressure of electrons which provides the main part of the total pressure. 
Furthermore, the neutrino emission associated with the electron capture reduces the 
thermal pressure of the core. Here the inner part of the core, which collapses nearly 
homologously with a subsonic infall velocity, constitutes the inner core. On the other 
hand, the outer region in which the infall velocity is supersonic constitutes the outer 
core. 

The collapse proceeds until the central part of the iron core reaches the nuclear 
density (~ 2 x 10 14 g/cm 3 ), and then, the inner core experiences the bounce. Because 
of its large inertia and large kinetic energy induced by the infall, the inner core 
overshoots its hypothetical equilibrium state. The stored internal energy of the 
inner core at the maximum compression is released through strong pressure waves 
generated inside the inner core. The pressure waves propagate from the center to the 
outer region until they reach the sonic point located at the edge of the inner core. 
Because the sound cones tilt inward beyond the sonic point, the pressure disturbance 
cannot propagate further and forms a shock just inside the sonic point. A shock wave 
is formed at the edge of the inner core and propagates outward. 

After this phase, the proto-neutron star experiences the quasi-static evolution 
phase. In this phase, the central value of density (the lapse function) increases 
(decreases) gradually, because the matter in the outer region falls into the proto- 
neutron star and because neutrinos are emitted carrying away the energy and lepton- 
number from the proto-neutron star. 

Figure [2] shows the radial profiles in the equator of the lepton fractions at the 
bounce. The central values of the electron, the electron-neutrinos, and the total- 
lepton fractions are « 0.32, 0.05, and 0.37, respectively. The electron-anti-neutrino 
fraction is almost zero through out the core because only very small amount of 
positrons exist due to the high degree of electron degeneracy 

6.3. Neutrino bursts and stall of shock 

As the shock wave propagates outward, the kinetic energy of the infall matter 
is converted into the thermal energy behind the shock. The conversion rate of infall 
kinetic energy may be estimated approximately as 



^heat ~ 47r^(p infall u 3 nfall /2) 
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Fig. 2. The radial profiles of the electron, v e -, i/ e -, and the total lepton fraction at the bounce. The 
results for the finer grid resolution (solid curve) and for the coarser grid resolution (the dotted 
curves) are shown together. The two results are almost identical. 
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~ 1.4 x 10 ergs/s : — ~ , (6-2) 

&/ \l00km) Vl0 9 g/cm 3 y V 0.2c / ' V ' 

where R s and p in f & \\ are radius of the shock wave and the density of infall matter, 
and we here recover the velocity of the light (c). Here, we assume that all the kinetic 
energy is converted to the thermal energy. 

At the same time, the shock wave suffers from the energy loss by the photo- 
dissociation of the iron to o-particles and free nucleons. The fraction of this energy 
loss is^ 

e diss ~ 1.5 x 10 51 ergs per O.1M . (6-3) 
Thus, the energy loss rate due to the photo-dissociation is 

Aiiss ~ M shock e diss ~ 1.1 x 10 53 ergs/s ( T7 ^—] ' ( Anfa " . ) ( ^) , (64) 
dlhh stock diss b/ ^ 100km y V.10 9 g/cm 3 y V 0.2c J ' v ; 

where M s hock ~ 47ri?g/9i n f a iit>i n f a n is mass-accretion rate to the shock front. 
The ratio of L heat to L diss is 

-^heat -i ~ / ^infall \ ^ / r , -\ 

z^" L Hw) • (6 - 5) 

Therefore the energy loss rate by the photo-dissociation will eventually overcome the 
hydrodynamic power, because the infall velocity, which is « {GM- lc / Rg) 1 ' 2 , decreases 
as the shock wave propagates outward. 
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Fig. 3. Time evolution of the neutrino luminosities. The results in the finer grid resolution (solid 
curves) and in the coarser grid resolution (dashed curves) are shown together. The two results 
are approximately identical until the convective phase sets in, whereas small disagreement is 
found in the convective phase. 

Furthermore, when the shock wave crosses the neutrino-sphere, spiky burst emis- 
sions of neutrinos, the so-called neutrino bursts, occur: Neutrinos in the hot post- 
shock region are copiously emitted without interacting core matter. Figure [3] shows 
the neutrino luminosity as a function of time calculated by^ 1 '^ 



The peak luminosity is L Ue ~ 4.5 x 10 53 ergs/s. This neutrino burst significantly 
reduces the thermal energy of the shock. Consequently, the shock wave stalls at 
~ 80 km soon after the neutrino burst. The peak luminosity and the shock-stall ra- 
dius agree approximately with the previous one-dimensional fully general relativistic 



When the shock wave stalls, negative gradients of the entropy per baryon and the 
total-lepton (electron) fraction appear because neutrinos carry away both the energy 
and the lepton number. Figure H] shows the radial profiles of the infall velocity, the 
density, the entropy per baryon, and the total lepton fraction in the equator. This 
figure clearly shows that negative gradients of the entropy per baryon and the total 
lepton fraction are formed above the neutrino sphere. As we shall see in § 16.41 such 
configurations are known to be unstable to convection, which is known to as the 
proto-neutron star convection. 
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Fig. 4. The radial profiles of the infall velocity, the density, the entropy per baryon, and the total 
lepton fraction at bounce, at 2 ms and 6ms after bounce. The results for the finer grid resolution 
(solid curves) and for the coarser grid resolution (the dotted curves) are shown together and 
they are shown to be approximately identical. 



6.4. Convective activities 

Let us investigate the stability of the envelope of the proto-neutron star following 
Lattimer and Mazurek!23 We consider the following parameter 



y2 _ 9cS 
P 



dp 
dr 



amb 



dp 
dr 



blob. 



(6-7) 



where g e Q is the effective gravitational acceleration defined to be positive in the 
negative radial direction, the subscript 'amb' refers to the ambient core structure, 
and 'blob' denotes the blob element which is under an isolated displacement. The 
condition iV 2 < implies that the structure is unstable to convective overturn (e.g. 
Ref.EZD). 

Assuming that the fluid elements maintain the pressure equilibrium with its 
surroundings, we have 
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(6-8 
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Fig. 5. Snapshots of the contours of the density (top left panels), the electron fraction Y e (top right 
panels), the entropy per baryon (bottom left panels), and the local neutrino energy emission 
rate (bottom right panels) in the x-z plane at selected time slices. 



Using this relation, Eq. (|6-7p is written as 
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Because the pressure is a function of the entropy per baryon, the density, and the 
lepton fraction, (dP/dr) am b is rewritten to 



^2 = 9eS I dp_ 



p \dP 



dP 

ds 



P,Yi 



*/amb 



dP 
dYi 



p.s 



dr J 



amb 



(6-10) 



Here, we also assume that the blob elements do not interact the ambient matters 
both thermally and chemically, i.e. ds = dY[ = for the blob. Then, we have 



d P / bi..b 



[dp 



(6-11) 



s,Y t 



Equation (|6-1U|) shows that when the pressure derivatives of given EOS ((dP/ds) p Y e 
and (dP/dYi) ps ) are positive, configurations with negative gradients of entropy and 
Y\ (N 2 < 0) are unstable. (Note that in the above treatment, we have ignored the 
dissociative effects caused by energy and lepton transports due to neutrinos.) Thus, 
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Fig. 6. Snapshots of the contours of gradients associated with the entropy per baryon 
(dP/ds)Y l , P (ds/dr) (right panels) and associated with the lepton fraction (dP/dYi) S:P (dYi/dr) 
(left panels) in the x-z plane at selected time slices. 



the negative gradients of the entropy per baryon and the total lepton fraction formed 
above the neutrino sphere lead to the convective overturn (the proto-neutron star 
convection). Indeed, convection occurs in our simulation. 

Figure [5] shows contours of the density, the electron fraction, the entropy per 
baryon, and the neutrino energy-emission rate. Convective motions are activated at 
about 8 ms after the bounce in the region located above the neutrino-sphere where 
the gradients of the entropy per baryon and Yi are imprinted (see Fig. 2]). At about 
10 ms after the bounce, the lepton rich, hot blobs rise to form 'fingers' (see in top 
left panel in Fig. Note that the neutrino energy emission rate in this finger is 
relatively higher than that in other region. This is responsible for the small hump 
seen in the time-evolution of neutrino luminosity (see Fig. [3|). Subsequently, the hot 
fingers expand to form 'mushroom structures', and push the surface of the stalled 
shock (see top right panel in Fig. [5]). At the same time, the lepton poor, colder 
matters sink down to the proto-neutron star (r < 20 km). The entropy per baryon 
just behind the shock increases to be s > Wks an d the stalled shock gradually moves 
outward to reach r ~ 200 km. As the hot, lepton rich matters are dug out from the 
region below the neutrino-sphere, the neutrino luminosity is enhanced (see Fig. 
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However, the energy released in the convective overturn is not sufficient to keep 
pushing the shock wave, and eventually, the shock stalls and turns to be a standing 
accretion shock (bottom two panels of Fig. [5]). All these features qualitatively 
agree with the previous multidimensional Newtonian simulations P^ 1 '® ESB ED 
A more detailed comparison with the previous simulations is given in § 16.51 

It will be interesting to investigate which gradient (entropy per baryon or elec- 
tron fraction) is more responsible for the convection. To see this, we calculate the 
gradients associated with the entropy per baryon (dP/ds)Y l , P (ds/dr) (right panels 
in Fig. [6]) and associated with the lepton fraction (dP/dYi) StP (dYi/dr) (left panels in 
Fig. [6]). This figure clearly shows that negative gradient of the entropy per baryon 
is more important for the convection activated promptly after the bounce. 

6.5. Comparison with the previous studies 

To check the validity of the code, the results presented in § 16. 2\ § 16.31 and § 16.41 
are compared with the previous simulations. 

6.5.1. Comparison of the results before the convection sets in 

We first compare our results with those in the state-of-the-art one-dimensional 
(ID) simulations in full general relativity,^'^'^'^ in which ID general relativistic 
Boltzmann equation is solved for neutrino transfer with relevant weak interaction 
processes. Because neutrino heating processes [y e + n — > p + e~ and v e + p — > 
n + e + ) are not included in the present implementation, and on the other hand, 
multidimensional effects such as convection cannot be followed in the one-dimensional 
reference simulations, we pay particular attention to comparing results during the 
collapse and the early phase (~ 10 ms) after the bounce (see results in § 16.21 and 

Our radial profiles of the lepton fractions at the bounce (see Fig. [2]) approxi- 
mately agree or at least are consistent with the previous simulations, implying that 
our code can correctly follow the collapse until the bounce. Also, the radial profiles of 
the infall velocity, the density, and the entropy per baryon just after the bounce show 
good agreements with the previous studies. No such good agreement was reported 
in the previous simulations^ '^D where simple leakage schemes based on the single 
neutrino-trapping density were adopted. Quantitatively, the negative gradients of 
the entropy per baryon and the lepton fraction are little bit steeper in the present 
simulation than those in ID full Boltzmann simulations. The reason may be partly 
because the transfers of lepton-number and energy are not fully solved in the present 
leakage scheme. Except for this small quantitative difference, the two results agree 
well. 

For validating a scheme for the neutrino cooling, agreement of the neutrino lumi- 
nosities with those by ID full Boltzmann simulation should be particularly checked 
because they depend on both implementations of weak interactions (especially elec- 
tron capture in the present case) and treatments of neutrino cooling (the detailed 
leakage scheme). Also, accurate computation of the neutrino luminosities is required 
for astrophysical applications, because neutrinos carry away the most of energy lib- 
erated during the collapse as the main cooling source and can be primary observable. 
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Our results, in particular the duration and the peak luminosity of the neutrino bursts, 
agree approximately with those in the previous simulations. Again, no such good 
agreement was reported in the previous simulation.®'^ 

The shock stall-radius is ~ 80 km. This value is consistent with (although 
slightly smaller than) that in Liebendorfer et al.® (-Rstaii ~ 85 km) and smaller 
than that in Sumiyoshi et al. 16) (i? s taii ~ 100 km). This is likely because in our 
leakage scheme, neutrino heating is not taken into account. 

To summarize, the results in the present simulation agree well with those in the 
previous ID Boltzmann simulations qualitatively. Quantitatively, the present results 
agree approximately with those in the previous ID Boltzmann simulations. We can 
obtain approximately correct results with a not computationally expensive scheme 
without solving the Boltzmann equation. Thus, the present code may be adopted, 
as a first step, to other multidimensional simulations such as the rotating stellar 
collapse to a black hole and mergers of compact binaries. 

6.5.2. Comparison of the results after the convection sets in 

In this section, we compare our results in the convective phase with those in the 
two-dimensional (2D) Newtonian s imulationP»™™®«™>E3 - m which 
a wide variety of approximations were adopted for the treatment of neutrinos. 

In the present simulation, we have found both the vigorous convective activities 
(the proto-neutron star convection) and the enhancement of neutrino luminosities 
due to the convection. These features agree approximately with those in the previous 
2D simulations with a fluid-like treatment of neutrinos^ 1 and with radial ray-by- 
ray, gray flux-limited diffusion approximation of neutrino transfers.® 4102141011 i n 
a spherically symmetric, gray flux-limited diffusion schemep^} D y contrast, only 
mildly active convection was found and no enhancement in the neutrino luminosities 
was observed. 

Note that the transport of energy and lepton number by neutrinos can flatten 
the negative gradients of entropy and lepton fraction, and as a result, the convec- 
tion will be suppressed. In purely hydrodynamic simulations without neutrino pro- 
cesses®''^^ 1 (using a postbounce core obtained in ID simulations with neutrinos), 
the proto-neutron star convection is strongly activated. In the radial ray-by-ray sim- 
ulations,® 4UE34Unj the transfer of neutrinos in the angular direction is not taken into 
account and the stabilizing effect is underestimated, resulting in the proto-neutron 
star convection with the enhancement of neutrino luminosities. In the spherically 
symmetric simulation,^--' the transfer of neutrino in the angular direction is as- 
sumed to occur fast enough to make the neutrino distribution function spherically 
symmetric, and consequently, the stabilizing effect is overestimated. 

Recently, Buras et al.^ 1 performed simulations with a modified ray-by-ray, 
multi-group scheme in which some part of the lateral components are included, and 
found that the proto-neutron star convection indeed sets in but has minor effects 
on the enhancement of the neutrino luminosities. Dessart et alP® performed simu- 
lations employing a 2D multi-group flux-limited diffusion scheme and found similar 
results as in Buras et al. Thus, although the proto-neutron star convection indeed 
occurs, its influence on enhancing the neutrino luminosities may be minor. The 
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Fig. 7. Gravitational wave quadrupole amplitude Ai due to the prompt convection as a function 
of post bounce time The results for the finer grid resolution (solid curve) and for the coarser 
grid resolution (the dotted curves) are shown together. 



strong convective activities and the enhancement of neutrino luminosities found in 
the present simulation should be considered as the maximum ones. 

Note that it is in intermediate regions (t„ ~ 1) that the stabilizing effect due 
to the neutrino transfer works efficiently: At higher density region with Tj, > 1, 
neutrinos cannot efficiently transport the energy and the lepton number due to the 
large opacities; At lower density region with t„ <C 1, on the other hand, neutrinos 
carry away the energy and the lepton number without interacting with the mat- 
ter. Therefore a careful and detailed treatment of the neutrino transfer is required 
to clarify the degree of the stabilizing effect and the convection, although such a 
computationally expensive sumulation is beyond the scope of this paper. 

The present result that the proto-neutron star convection occurs qualitatively 
agrees with the recent simulations with detailed neutrino transfer liQSI JZ2J jf simula- 
tions are perfomed keeping in mind that the stabilizing effect due to the neutrino 
transfer is not taken into account in the present scheme, the present code will be 
acceptable to explore the the rotating stellar collapse to a black hole and mergers of 
compact binaries. 

6.6. Gravitational radiation 

Associated with the convective motions, gravitational waves are emitted. The 
gravitational waveforms are computed using a quadrupole formulaP^ In quadrupole 
formulae, only the +-mode of gravitational waves with I = 2 and m = is nonzero 
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Fig. 8. The frequency spectra of the characteristic gravitational-wave strain due to the prompt 
convection. The results for the finer grid resolution (solid curve) and for the coarser grid 
resolution (the dotted curves) are shown together. 



in axisymmetric spacetime and is written as 

7 quad Izz\trct) ~ Ixx{tret) ■ 2 a — -^2^) ■ 2 a frin\ 

/il = sin u = sin 0, (6-12) 

r r 

where ly denotes a quadrupole moment, ly its second time derivative, and t ret a 
retarded time. In fully general relativistic and dynamical spacetime, there is no 
unique definition for the quadrupole moment and nor is for Ijj. Following Shibata 
and Sekiguchi we choose the simplest definition of the form 



J y = J p.jVd'j. (6-13) 
Then, using the continuity equation, the first time derivative can be written as 

I ij = J ft ( V + x i v j )d 3 x, (6-14) 

and Iij is computed by the finite differencing of the numerical result for Ijj. In the 
following, we present A2, which provides the amplitude of a given mode measured 
by an observer located in the most optimistic direction (in the equatorial plane) . We 
also calculate the characteristic gravitational-wave strain 



w/) = V^§z^' (645) 
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Fig. 9. Snapshots of the contours of \/—N 2 /2tv, in the x-z plane at selected time slices. 



where 



dE 



Svr 2 c 3 _ /2 



Mf) 



15 G' 

is the energy power spectra of the gravitational radiation and 



Mf) = J A 2 (t)e 2 ^ t dt. 



(6-16) 



(6-17) 



Figure [7] shows Azii). Because the system is initially spherically symmetric, 
no gravitational radiation is emitted before the onset of the convection. When the 
proto-neutron star convection sets in at ~ 10 ms after the bounce, gravitational 
waves start to be emitted. The peak amplitudes are A 2 ~ 100 cm. After the peak is 
reached, gravitational waves generated by the smaller-scale convective motions are 
emitted with A 2 ~ 50 cm. 

Figure [8] shows the spectra of /i c har due to the convective motions. In contrast 
to the spectra due to the core bounce (e.g. Refs. I20|) and l30j) ). there is no dominant 
peak frequency in the power spectra. Instead, several maxima for the frequency 
range 100-1000 Hz are present. Note that for gravitational waves due to the core 
bounce, the characteristic peak frequency is associated with the bounce timescale 
of the core. The effective amplitude of gravitational waves observed in the most 
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optimistic direction is /i c har ~ 6-8 x 10 for an event at a distance of 10 kpc, which 
is as large as that emitted at the bounce of rotating core collapsed 

To check that gravitational waves are indeed originated by the convective mo- 
tions, we calculate the frequency \/—N 2 /2ir (see Eq. (|6-7p ) as shown in Fig. [9l This 
frequency is in good agreement with the gravitational-wave frequency implying that 
gravitational waves are indeed due to the convective activities. 

Miiller and JanksP^ 1 investigated gravitational waves due to the convective mo- 
tion inside the proto-neutron star. It is interesting to compare our results with theirs. 
They adopted a post-bounce model of Hillebrandt^^* They put an inner boundary 
at radius r- m = 15 km and assumed the hydrostatic equilibrium there. They do not 
include neutrino transfer while a sophisticated EOS is adopted. They found quali- 
tatively similar results to ours. According to their results, the maximum amplitude 
of the quadrupole mode is A<i ~ 100 cm, which agrees well with our results. The 
spectrum of the gravitational-wave strain has several maxima for / = 50-500 Hz 
with the maximum value of /i c har ~ 3 x 10~ 21 . The peaks in /i c har are distributed for 
higher frequency side in our results probably due to the general relativistic effects. 
We note that a similar general relativistic effect is observed for gravitational waves 
at the bounce phaseP^ These facts show that for deriving quantitatively correct 
spectra of gravitational waves, fully general relativistic simulations are necessary. 

6.7. Numerical accuracy 

In Figs. [TH3]we show the results both in the higher resolution (solid curves) and 
in the lower resolution (dashed curves). The radial profiles of the two resolutions 
are almost identical, showing that convergent results are obtained in the present 
simulation (see Fig. In the time evolution of neutrino luminosities (see Fig. [3j), 
the two results are almost identical before the convective activities set in. In the 
later phase, on the other hand, the two results show slight disagreement. Because 
the convection and the turbulence can occur in an infinitesimal scale length, the 
smaller-scale convection and turbulence are captured in the finer grid resolution. 
However, the influence of the grid resolution on the neutrino luminosities is minor 
because the convection and turbulence are strongly activated in the region above the 
neutrino sphere (see the contours of the electron fraction and the entropy in Fig. 
[5]). On the other hand, most of the neutrinos are emitted from the region inside the 
neutrino sphere (see the contour of the local neutrino energy emission rate in Fig. 

The effect of the grid resolution can be seen in gravitational waves. In Figs. [7] 
and[S]we show the quadrupole mode Azif) and the characteristic strain /i c har(/) both 
in the higher resolution (solid curves) and in the lower resolution (dashed curves). 
After the formation of the lepton-rich, hot finger at ~ 10 ms after the bounce (see 
§ 16. 4p . convective activities set in. Then, disagreement of A<i(t) between the finer 
and the coarser grid resolutions becomes noticeable (see Fig. [7J. The characteristic 
peaks of /i c har(/) i n higher frequencies (/ ~ 200-500 Hz) are more prominent (see 
Fig. [HJ). It is likely to be because the smaller-scale turbulant motions are captured 
in the finer grid resolution. 

To check the accuracy of our numerical results, the violation of the Hamiltonian 



Full GR simulation with microphysics 



33 




Fig. 10. Evolution of the averaged violation of the Hamiltonian constraint (upper panel) and baryon 
mass conservation (lower panel). 



constraint are calculated, which is written as 
H = -8ip~ 5 



(6-18) 



where tp = an d A denotes the Laplacian with respect to jij. In this paper, the 
averaged violation is defined according to® 1 



ERROR = — / p*\V\d 3 x, 
M* J 

where M* is the rest-mass density of the core (see Eq. 



(6-19) 
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(6-20) 



Namely, we use as a weight factor for the average. This weight factor is 
introduced to monitor whether the main bodies of the system (proto-neutron stars 
and inner cores), in which we are interested, are accurately computed or not. 

We display the time evolution of the Hamiltonian-constraint violation and the 
conservation of the baryon mass of the system in Fig. [TUl Several discontinuous 
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changes in the Hamiltonian-constraint violation and the conservation of the baryon 
mass originate from the regridding procedures in which some matters of the outer 
region are discarded. 

Before the bounce, the baryon mass is well conserved and the Hamiltonian- 
constraint violation is very small as ~ 10" 4 . After the bounce, the violation of 
the baryon-mass-conservation and the Hamiltonian constraint is enhanced due to 
the existence of shock waves where the hydrodynamic scheme becomes essentially 
a first-order scheme. The convergence of the baryon-mass-conservation and the 
Hamiltonian-constraint violation also becomes worse in the convective phase. How- 
ever, the degree of violation of the Hamiltonian constraint and the baryon-mass- 
conservation is small and we may believe that the numerical results obtained in the 
paper are reliable. 

§7. Summary and Discussion 

7.1. Summary 

In this paper, we present a fully general relativistic hydrodynamic code in which 
a finite-temperature EOS and neutrino cooling are implemented for the first time. 
Because the characteristic timescale of weak interaction processes t wp ~ | Y e /Y e | (WP 
timescale) is much shorter than the dynamical timescale t^ yn in hot dense matters, 
stiff source terms appear in the equations. In general, an implicit scheme may be 
required to solve themJES However, it is not clear whether implicit schemes do work 
or not in the relativistic framework. The Lorentz factor is coupled with the rest- 
mass density and the energy density. The specific enthalpy is also coupled with the 
momentum. Due to these couplings, it is not straightforward to recover the primitive 
variables and the Lorentz factor from conserved quantities. Taking account of these 
facts, we proposed an explicit method to solve all the equations noting that the 
characteristic timescale of neutrino leakage from the system ti ea k is much longer 
than t wp and is comparable to tdyn- 

By decomposing the energy-momentum tensor of neutrinos into the trapped- 
neutrino and the streaming-neutrino parts, the hydrodynamic equations can be 
rewritten so that the source terms are characterized by the leakage timescale ti ea k 
(see Eqs. (|3-8p and (j3-5|) V The lepton- number conservation equations, on the other 
hand, include the source terms characterized by the WP timescale. Taking account 
of these facts, limiters for the stiff source terms are introduced to solve the lepton- 
number conservation equations explicitly (see § 13. 2|) . In the numerical relativistic 
hydrodynamics, it is required to calculate the primitive variables and the Lorentz 
factor from the conserved quantities. In this paper, we develop a robust and stable 
procedure for it (§ 15. 3p . 

To check the validity of the numerical code, we performed a simulation of spheri- 
cal stellar core collapse. As initial conditions, we adopted the 15Mq spherical model 
with the solar metallicity computed by Woosley et alJ® After the shock formation 
and propagation, the shock wave suffers from the severe reduction of its energy due 
to neutrino burst emission when the shock wave passes the neutrino-sphere. Even- 
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tually, the shock wave stalls soon after it passes through the neutrino sphere. The 
neutrino burst makes negative gradients of the entropy and Y[ above the neutrino 
sphere. Because such configuration is convectively unstable, vigorous convective mo- 
tions are induced. All these properties agree qualitatively with those by the resent 
2D Newtonian simulations.'^' 1 ''^' 

We also compared our results with those in the previous simulations. Before 
the convection sets in, we compare our result with those in the state-of-the-art ID 
Boltzmann simulations in full general relativity QSJJSDJISHniJ As shown in this paper, 
the radial structure of the core and the neutrino luminosities agree qualitatively well 
with those in their simulations. Quantitatively, they also agree approximately with 
the previous results. 

After the convection sets in, we compare our result with those in 2D Newtonian 
simulations .128JI ,[99ji ,liooj» ,lioij ,[23J) ,H02j) ,ii03ji ,l72ji Q ur resu ^ that the proto-neutron star con- 
vection occurs agree qualitatively with that in the previous simulations ^ ^Sb .LLQ3J ,E3 
However, quantitative properties show disagreement because the transfer of neutri- 
nos are not fully solved in the present scheme. Note that the transport of energy and 
lepton-number by neutrinos can flatten the negative gradients of entropy and lepton 
fraction, stabilizing the convection. Therefore the convective activities obtained in 
the present simulation should be considered as the maximum ones. 

If we keep in mind the above facts and note the good agreements of the radial 
structure and neutrino luminosities, the present implementation will be applied to 
simulations of rotating core collapse to a black hole and mergers of binary neutron 
stars as a first step towards more sophisticated models. A detailed treatment of the 
neutrino transfer is required to determine the degree of stabilizing effect, but this is 
far beyond the scope of this paper. 

Gravitational waves emitted by the convective motions are also calculated. The 
gravitational-wave amplitude is ~ 3 x W~ 21 for an event of the distance 10 kpc. Re- 
flecting the contributions of multi-scale eddies with characteristic overturn timescale 
1-10 ms, the energy power spectrum shows several maxima distributed in / ~ 100- 
1000 Hz. We compare our results with those in Miiller and JankcP^ 1 in which a similar 
calculation (but in Newtonian gravity) is performed. The maximum amplitude of 
gravitational waves in our results agrees well with that in Miiller and Janka. The 
several maxima in the energy power spectrum are distributed at higher-frequency 
side in our results due to the general relativistic effects, showing that fully general rel- 
ativistic simulations are necessary for the accurate calculation of gravitational-wave 
spectra. 

7.2. Discussions 

Because the present implementation of the microphysics is simple and explicit, 
it has advantage that the individual microphysical processes can be easily improved 
and sophisticated. For example, the neutrino emission via the electron capture 
can be easily sophisticated as follows. To precisely calculate the electron capture 
rate, the complete information of the parent and daughter nuclei is required. In 
EOSs currently available, however, a representative single-nucleus average for the 
true ensemble of heavy nuclei is adopted. The representative is usually the most 
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abundant nucleus. The problem in evaluating the capture rate is that the nuclei 
which cause the largest changes in Y e are neither the most abundant nuclei nor the 
nuclei with the largest rates, but the combination of the two. In fact, the most 
abundant nuclei tend to have small rates because they are more stable than others, 
and the fraction of the most reactive nuclei tend to be smallP^^ Assuming that 
the nuclear statistical equilibrium (NSE) is achieved, the electron capture rates under 
the NSE ensemble of heavy nuclei may be calculated for a given set of (p, Y e , T). 
Such a numerical rate table can be easily employed in the present implementation. 

Also, the neutrino cross sections can be improved. As summarized in Ref. 1108]) . 
there are a lot of corrections to the neutrino opacities. Note that small changes 
in the opacities may result in much larger changes in the neutrino luminosities, 
because the neutrino energy emission rates depend strongly on the temperature, and 
the temperature at the last scattering surface (t„ ~ aT 2 ~ 1) changes as T ~ c -1 / 2 . 
Although the correction terms are in general very complicated, it is straightforward 
to include the corrections into our code. Note that the corrections become more 
important for higher neutrino energies. Therefore, the correction terms might play a 
crucial role in the collapse of population III stellar core and the formation of a black 
hole, in which very high temperatures (T > 100 MeV) will be achieved. A study 
to explore the importance of these corrections in the case of black hole formation is 
ongoing. 

As briefly described in the introduction, one of the main drawbacks in the present 
implementation of the neutrino cooling is that the transfer of neutrinos are not 
solved. Although fully solving the transfer equations of neutrinos is far beyond the 
scope of this paper, there are a lot of rooms for improvements in the treatment of 
the neutrino cooling. For example, the relativistic moment formalism '^323 [ n 
particular the so-called Ml closure formalism, may be adopted. For this purpose, a 
more sophisticated treatment of the closure relation for P a a is required. We plan to 
implement a relativistic Ml closure formalism for the neutrino transfer in the near 
future. 

To conclude, the present implementation of microphysics in fully general rela- 
tivistic, multidimensional code works well and has a wide variety of applications. We 
are now in the standpoint where simulations of stellar core collapse to a black hole 
and merger of compact stellar binaries can be performed including microphysical 
processes. Fruitful scientific results will be reported in the near future. 
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Appendix A 

Electron and positron captures 

In this section, we briefly summarize our treatments of electron and positron 
captures which are based on Ref. 66) and give the explicit forms of 7^, j$ e , Q1 c e , 
and Qpl appeared in Eqs. (|3-27j) . (|3-28p . and (|3-30p . for the purpose of convenience. 

A.l. The electron and positron capture rates 7^ and 7?g 

The 'net' electron fraction is written as Y e = Y_ — Y + where Y_ (1+) denotes 
the number of electrons (positrons) per baryon including pair electrons. Then the 
electron-neutrino number emission rate by the electron capture and the electron- 
anti-neutrino number emission rate by the positron capture are given by 

7 local = _y_ = _ { yf + yhj ^ ^ 

^ = -Y + = -(Yf + Y<t), (A-2) 

where the electron and positron capture rates are decomposed into two parts, capture 
on by free nucleons (with superscript /) and on heavy nuclei (with superscript h). 
In the following, we present the explicit forms of y£, Y+, Yh, and Y+. 

A. 2. Capture on free nucleons Y* 

The electron capture rate (including the contribution of the inverse reaction of 
the neutrino capture) on free nucleons (y£) is given by 

yL = X n \ UeCj - X p \ ec ' f , (A-3) 

where A ec '^ is the specific electron capture rate on free protons, X VeC ^ is the specific 
electron-neutrino capture rate on free neutrons, and X p and X n are the mass fraction 
of free proton and neutron, respectively. Based on a balance argument,^ 1 one can 
show that \ VeC ^ is related to A 60 '-' by 

\ VeC, f = exp L ve - Ve - ^\ \ ec >f, (A-4) 

where r\ ve and r\ e are the chemical potentials of electron neutrinos and electrons in 
units of ksT and 5m = (m n — m p )c 2 . Furthermore, we use the following relation for 
non-degenerate free nucleons, 

where r] n and r] p are the chemical potentials of free neutrons and free protons in units 
of ksT. Then we obtain 

Yl = [exp (rj ve - r/ e + r/ n - r, p ) - 1] X p X ec ' f . (A-6) 

The positron capture rate (including the contribution of the inverse reaction) on 
free nucleons is similarly given by 

Yl = X p X v ^f - X n X^f = [exp ( VPe + Ve + rj p - Vn ) - 1] X n X^ , (A-7) 
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where r]p e is the chemical potential of electron-anti- neutrinos in units of kpT, A pc is 
the specific positron capture rate on free neutrons, and \ UeC ^ is the specific electron- 
anti- neutrino capture rate on free protons. 

A. 3. Capture on heavy nuclei Y h 

The electron capture rate (including the contribution of the inverse reaction of 
the neutrino capture) on a heavy nucleus of mass number A (Yh) is given by^ 



yn 



X P X ec,h 



(A- 



where \ ec > h is the specific electron capture rate on the parent nucleus (mass fraction 
Xp), \ UeC ' h is the specific electron neutrino capture rate on the daughter nucleus 
(mass fraction Xp>), and A is the atomic mass of the parent and daughter nuclei. 
In the present simulations, we set Xp> = Xp = Xa- Then, under the assumption 
of nuclear statistical equilibrium, one may approximate the capture rate on heavy 
nuclei as J® 

Y h « [exp ( Vve -Ve + Vn- V P ) - 1] ^A 6 ^. (A-9) 



A 



Similarly, the positron capture rate (including the contribution of the inverse 
reaction) on heavy nuclei (Yj£) is given by 



Yl 



~A X 



Xp 
A 



X pc ' h « [exp (i] Pe + n e + T] p - r} v 



1] ^\ pc > h . 
1 A 



(A40) 



A. 4. The specific capture rate A 

The specific electron and positron capture rates on free nucleons and on heavy 
nuclei and are written in the same form as®! 



A' 



ec,/ 



A' 



ec,h 



In 2 
In 2 

cm? 



ec,/ 



ec,h 



X pcj 



In 2 



In 2 



I pcJ , 

jpc,h 



(A-ll) 
(A-12) 



where I ec, $ and I pc ^ are the phase space factors for the electron and positron cap- 
tures on free electrons, and I ec,h and I pc > h are those on heavy nuclei. (/i) e ff's are 
the effective ft- values introduced by Fuller et al.j® which is essentially the same as 
the square of the nuclear transition matrix. 
The phase space factors are given by 



jec,/ . 
jpc,/ . 
rec,h . 
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1 
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air], (A-14) 
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dr], (A-16) 



1 + e »?-^e+C pc ' h 



where £ ec '^ ; C pc '^> C ec,h , and C pc ' ft ' are the nuclear mass-energy differences for electron 
capture and positron capture in units of ksT. The nuclear mass-energy differences 
for capture on free nuclei are given by 

r ,f = _ C pc,f^ rjp _ rjn (A . 17) 

We follow Fuller et alP^ for the nuclear mass-energy differences for capture on heavy 
nuclei: In the case of N < 40 or Z > 20 (referred to as 'unblocked' case), we set 

C^ h = -C pc ' h ^V P -Vn. (A-18) 

In the case of N > 40 or Z < 20 (referred to as 'blocked' case), on the other hand, 
we set 

-» + * + (A-20) 

Then, the threshold value of the electron and positron captures is given by rjo = 
m e c 2 /{k,BT) for £ > —m e c 2 /(kBT) and tjq = |£| for ( < —m e c 2 /(kBT) where we 
have dropped the superscripts 'ec', 'pc', '/', and '/i' in £ for simplicity. 

The effective /t-value of electron or positron capture on free nuclei is given by 
(e.g. Ref.EHD 

Icgio^S' = logio(/i) e P ff f ~ 3.035. (A-21) 

We follow Fuller et alP^ for the effective /t-value of capture on heavy nuclei, who 
proposed to use 

{3.2 unblocked r/ e < \( ec > h \ 
2.6 unblocked Ve > \( ec > h \ , (A-22) 
2.6 + fr 9 - blocked 

r 3.2 unblocked r/ e < |C pc ' h | 
logio(/*)eff H ~ < 2 - 6 unblocked r? e > |C pc ' ft | , (A-23) 
[ 2.6 + frp blocked 

where T 9 = T/(10 9 A). In this expression, the thermal unblocking effectP^ 1 is readily 
taken into account. In the thermal unblocking, it costs ~ 5.13 MeV to pull a neutron 
out of a filled orbital I/5/2 and place it in the gd-shell!^ 

A. 5. Energy emission rates Q^ c e and Q^ c e 

The neutrino energy emission rates associated with electron and positron cap- 
tures in units of m e c 2 s _1 are given by^ 

T6C TV)C 

7r ec = ln2— — , vr pc = ln2— w , (A-24) 

</*>s </*>s 
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where the phase space factors are given by 



(k B T\ l 
\m e c 2 J 



jpc _ (hsLx 



\m e c 2 J 



'Id 

DO 



> I" 



1 



1 + e 5 ?"^ 
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1 + e V+Ve 



1 



1 _ e rj-Vve+< c 

1 



1 _ e V-^e+( PC 



dr], (A-25) 
drj. (A-26) 



In Eqs. (|A-24p - (|A-26p . we have dropped the superscripts '/' and , h' in 7r ec , 7r pc , J ec , 
J pc , (ft)th C ec , and C pc for simplicity. 

The average energy of electron neutrinos produced by electron and positron 
captures is defined, in units of m e c 2 , as 



rec Tec 
If \ ec — l f - \P C — _ 



(A-27) 



Then, the local neutrino energy emission rates by the electron and positron captures 
per unit volume is given by 
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ec 
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X p (e ve )^\ c ^ + ^A{f ve )^ h \^ h 

Appendix B 

- Neutrino pair processes 



(A-28) 
(A-29) 



In this section, we briefly summarize our treatment of pair processes of neutrino 
emission and give the explicit forms of 7^, 7^, 7? e p e mS > 7^X' ^JT ' Qlj? e > 

Ql% Q„ r J? S ' QltL Qtl, and Qu'JT appeared in Eqs. & and 

(|3-30p . for the purpose of convenience. 

B.l. Electron-positron pair annihilation 

We follow Cooperstein et al.^J for the rate of pair creation of neutrinos by the 
electron-positron pair annihilation. The number emission rate of v e or v e by the 
electron-positron pair annihilation can be written as 



Pair 

WeVe 



m 



u Cvetl o-qc {k B Tf 



p 367T 4 m 2 c 4 (he) 



F 3 ( Ve )F 3 (-rj e ){Uo C k)l% 



(B-1) 



where a » 1.705 x 10- 44 cm- 2 and C* % = (C v - C A ) 2 + (Cy + C A ) 2 with C v = 
g + 2sin 2 #v^ and C A = 5- The Weinberg angle is given by sin 2 9\y ~ 0.23. Using 
the average energy of neutrinos produced by the pair annihilation, 



/ \pair _ f F 4 (r/ e ) F 4 (-77e) 

{ve9e) ~ 2 [fM + F 3 (-rj e ) 



(B-2) 
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the blocking factor (block) £ a p is evaluated as 



(block) 



pair 
ve 



1 + exp T], 



pair 



k B T 



1 + exp rj ue 



pair 



-1 -1 



The associated neutrino energy emission rate by the pair annihilation is given by 

(B-4) 



n pair _ _^_„.pair / _ vpair 



m,, 



Similarly, the number emission rate of ^ or by the electron-positron pair 
annihilation and the associated energy emission rate are given by 



.-~<pair /, rp\g 
pair m u ^v x v x °~0 C \ K B-L ) „ / \ t-, / \ /1 , i \pair 
— ^3(?7e)-F3(-?7e)(block)^ 



'fx 



p 367T 4 m|c 4 (/ic) 

^pair _ P pair / _ vpair 
^iu x u x lv x v x \ b VxVxl 1 



(B-5) 
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where C^p^ = (Cy — Ca) 2 + (Cy + Ca — 2) 2 . The average neutrino energy and the 
blocking factor are given by 



and 

(block)Pf 



1 + exp i] ux 



1 _ \pair _ I _ \pair 

I _ \pair 
V-VxVxl 



(B-7) 



1 + exp 77, 



( e v x v x ) 



pair 



(B-8) 



where note that r\y X = r/ ux . 

B.2. Plasmon decay 

We follow Ruffert et al. — ' for the rate of pair creation of neutrinos by the decay 
of transversal plasmons. The number emission rate of v e or v e can be written as 



pias = run C 2 V a c {k B Tf 6 +7 „)( b lock) p 



plas 



(B- 



where a^ ne « 1/137 is the fine-structure constant and 7 P rj 2Y / (ctfi ne /97r)(7r 2 +3rj e ). 
The blocking factor is approximately given by 

~! r / ' \plas \ "i ~1 



(block) |J? 
where 



1 + exp r) v 



plas 



k R T 



1 + exp 77, 



k R T 



(ei/ e i/ e ) 



plas 



2 + 



l + l7p 



(B40) 



(B41) 



is the average energy of neutrinos produced by the plasmon decay. The associated 
neutrino energy emission rate is given by 



n plas _ P _plas / _ vplas 



(B-12) 
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Similarly, the number emission rate of v x or v x by the plasmon decay and the 
associated energy emission rate are given by 



plas 



p 1927r 3 afi ne m 2 c A (he 

f)plas P plas / _ vplas 



TO,, 



(B44) 



where the average neutrino energy is (e^p^) 15 = (e I/e p e ) plas and the blocking factor 
is given by 
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B.3. Nucleon-nucleon bremsstrahlung 

We follow Burrows et al.^ for the rate of pair creation of neutrinos by the 
nucleon-nucleon bremsstrahlung radiation. They derived the neutrino energy emis- 
sion rate associated with the pair creation of v x or v x by the nucleon-nucleon 
bremsstrahlung radiation without the blocking factor: 

28 
~3 



Q^X S '° = 3.62 x 10°C 



' I Xf, ■ A'; • ==-X n X p )ft 



(k B T\^ 

1 1 K^VxVx) 



\m e c 2 J 



where £ Brems ~ 0.5 is a correction factor and the average energy is 
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(B-16) 
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We multiply Q^ V p ms, ° by the blocking factor, 
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to obtain the 'blocked' neutrino energy emission rate 

<er = QSr'°(biock)^ s - 

The number emission rate of u x or v x is readily given by 
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(B-20) 

Noting that the weak interaction coefficients of the bremsstrahlung radiation 
ar giiiji (1 — Cy) 2 + (1 — Ca) 2 for the pair creation of v x v x and C y + C\ for the pair 
creation of v e v e , the number emission rate and the associated energy emission rate 
for v e or v e may be written as 
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Appendix C 

Neutrino diffusion rates 

We follow Ref. [62]) for the diffusive neutrino-number emission rate j dlS and the 
associated energy emission rate Q dlS appeared in Eqs. (|3-25p and ()3-26p . and present 
the explicit forms of them in the following for convenience. An alternative definition 
of the diffusion rates are found in Ref. I63|) . 

C.l. Neutrino diffusion rates 

To calculate the neutrino diffusion rates 7^ lff and Q dlS , we first define neutrino 
diffusion time. In this paper, we consider cross sections for scattering on nuclei 
(cr^), on free protons (c^,), and on free neutrons (c^), and for absorption on free 
nucleons for electron neutrinos and a^Z for electron anti-neutrinos. 

Ignoring the higher order correction terms in neutrino energy E u , these neutrino 
cross sections can be written in general as 

a(E u ) = E 2 v a, (C-l) 

where a is a 'cross section' in which E% dependence is factored out. In practice, the 
cross sections contain the correction terms which cannot be expressed in the form 
of Eq. (|C-lj) . We take account of these correction terms, approximating neutrino- 
energy dependence by temperature dependence according to 

The opacity is written as 

k{E v ) = Y, K i( E u) = E 2 V S £ J ~^ = E% (C-3) 
and optical depth is calculated by 

T (E V ) = J K{E u )ds = E l J R ds = Elf. (C-4) 
Then, we define neutrino diffusion time by 

T^(E U ) ee a diS ^^r(E u ) = E^t = E^ , (C-5) 

C CK 

where the distance parameter Ax{E v ) is given by 

Ax{E v ) = ^1 (C-6) 

Note that f dm can be calculated only using matter quantities. a diff is a parameter 
which controls how much neutrinos diffuse outward and we chose it to be 3 following 
Ref. [63]) . For a larger value of a dlS , corresponding neutrino emission rate due to 
diffusion becomes smaller. 
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Finally, we define the neutrino diffusion rates by 

diff _ m u f n u (E u ) 1 m u 4ncg u k 



-<° s T 1 TM) dE " = ^TW TF ° M - (C ' 7) 

diS _ f E u n u (E u ) 1 Ancg v k 2 

Qv = J T^(E V ) dK ~ ^Jhc?-* T Fl{jlv) - (C ' 8) 

C.2. Summary of cross sections 

In this subsection, we briefly summarize the cross sections adopted in the present 
neutrino leakage scheme. 

C.2.1. Neutrino nucleon scattering 

The total u-p scattering cross section a p for all neutrino species is given by 

< P = f [i°y ~ !) 2 + ^\{Ca - 1) 2 ] , (C-9) 

where is the axial-vector coupling constant qa ~ —1.26. On the other hand, the 
total v — n scattering cross section a n is 

2 



sc = ^0 ( _Ev_ 

un 16 V mlc 2 



C.2.2. Coherent scattering of neutrinos on nuclei 

The differential cross section for the v-A neutral current scattering is written 



as^SJ 



da B J ao ( E v \ 2 2 h - 2 



dn-^\m e c*> » L [>VCV l+ C LO s]-(5, u ),l + cos„. ,CM1> 

where is the azimuthal angle of the scattering and 

2Z 

W = l-— (l-2sin 2 w )- (C-12) 

./I 

(<5ion)) Clqsj and Cfrp are correction factors d ue to the Coulomb interaction between 
the nucleip^ 1 due to the electron polarizationpSJ anc j d ue i Q the finite size of heavy 
nucleiP3) Because it is known that the correction factor Clos is important only for 
low neutrino energies,^ we consider only (<Si on ) and Cff- 

The correction factor due to the Coulomb interaction between the nuclei is given 

by 

(S ion ) = 7 / dcos#(l + cos0)(l - cos0)5 ion . (C-13) 

Itoh et alP^ presented a detailed fitting formula for the correction factor. However, 
the fitting formula is very complicated, we use a simple approximation based on 
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Ref. rTTTj) . 

<Si on can be written approximately as^^ 
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Fig. 11. Comparison of the ion-ion correction term, as a function of x = E v ai /(he), between our 
simplified estimation and a detailed fitting formula by Itoh et al. From the top to the bottom, 
the curves are for f = 1, 2, 5, 10, 20, 40, 80, 125, 160. 



for (qaj 1), where q = (2E U /Hc) sin(0/2), aj = (4^71^/3) ~ 1//3 is the ion-sphere 
radius, ua is the number density of a nucleus, and r = (Ze) 2 /(aiksT) is the con- 
ventional parameter that characterizes the strongness of the Coulomb interaction. 
f(r) is given by^ 

0.3414ir 1/4 + 0.05484r~ 1/4 . (C-15) 



f(r) w 0.73317 - 0.39890r 
The integration approximately gives for x 
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(C-16) 

To use this expression for the case of x > 1, we set the maximum value as (5; on ) = 
max(l, (<Sion)) where (<Si on ) = 1 corresponds to the case without the correction. Note 
that x can be larger than unitj0- In this case, the simple approximation on average 
underestimates the effect of the Coulomb interaction between the nuclei (see Fig. 
□ID- 
Figure [11] compares the correction term as a function of the parameter x calcu- 
lated by our simple approximation and by a detailed fitting formula by Itoh et alP^ 1 
For smaller values of r, the disagreements become prominent. Note that the typical 
value of r inside the collapsing core with T ~ 1 MeV, p ~ 10 12 g/cm 3 , A = 56 and 
Z = 26 ( 56 Fe) is T ~ 47. 



x ~ (_E„/27MeV)(pi2/j4) 1//3 where pi2 is the rest-mass density in units of 10 12 g/cm 3 
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C.2.3. Absorption on free neutrons 

The total cross section of the absorption of electron neutrinos on free neutrons 
is given by^^ 

2 



cr 
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E v + A 



up 
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E v + A 



up 
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(C-17) 



where Z\ np = m n c 2 —m p c 2 , and Wm is the correction for weak magnetism and recoil 
which is approximately given by 



W M ~ 1 + 1.1- 



m n c 
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(C-18) 



Similarly, the total cross section of the absorption of electron anti-neutrinos on 
free protons is given by^ v 



ab 
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l + Sg 2 A \ (E D -A 
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where W M is approximately given by 
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